#include <math.h>
#include "eno.h"
#include "fluids.h"

#include "eno.xtn"
#include "fluids.xtn"

extern void handle_triple();
extern void invert_matrix();
extern double density_above,density_below,density_middle;
extern double compute_curvature(),compute_delta(),
              compute_heaviside(),compute_surface_force();

static logical use_corner_difference=TRUE;

double get_contribution(in_distance,cutoff)
double in_distance,cutoff;
{
  if (in_distance <= -cutoff)
   return 1.0;
  else if (in_distance >= cutoff)
   return 0.0;
  else
   return 0.5*(1.0-sin(pi*in_distance/(2.0*cutoff)) );
}    /* get_contribution */

/* compute smoothed value given value of 3 level set functions.
    input : in_distance,distance2,distance3 - level set functions
            value_above_value_below,value_middle - some physical value of the fluids
             (e.g. viscosity or density)
            value_spread - amount to smooth
*/
double smooth_front(in_distance,distance2,distance3,
                 value_above,value_below,value_middle,
                 value_spread)
double in_distance,distance2,distance3,
       value_above,value_below,value_middle,value_spread;
{
double cutoff,d1s,d2s,d3s;

  cutoff = value_spread*h;
  if (! use_bmo)
   handle_triple(&in_distance,&distance2,&distance3,cutoff);
 
  d1s=get_contribution(in_distance,cutoff);

  if ((three_phase())&&(feed_back()) )
   {
    d2s=get_contribution(distance2,cutoff);
    d3s=get_contribution(distance3,cutoff);
    /*
    printf("d1s,d2s,d3s,vb,va,vm%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n",
      d1s,d2s,d3s,value_below,value_above,value_middle);
    */
    return (d1s*value_below+d2s*value_above+d3s*value_middle)/
           (d1s+d2s+d3s);
   }
  else
   {    /* code if not three phase */ 
    d2s=1.0-d1s;
    return (d1s*value_below+d2s*value_above);
   }   /* code if not three_phase */
}    /* smooth_front */

/* compute the density given the value of the 3 level set functions */
double cvt_density(in_distance,distance2,distance3)
double in_distance,distance2,distance3;
{
  return smooth_front(in_distance,distance2,distance3,
            density_above,density_below,density_middle,
            density_spread); 
}

/* surface tension update for pressure gradient update 
   only works for 2-fluid version */
void tensionforce(ix,iy)
int ix,iy;
{
double curv,denval,delta_function,holdval;
double distance_dd,dx,dy,out_dist;
int dnum;

  dnum=1;
  delta_function=compute_delta(gdt(dnum,ix,iy));
  curv=compute_curvature(ix,iy,&distance_dd,&dx,&dy,
   &out_dist,dnum); 
  denval=cvt_density(gd1(ix,iy),gd2(ix,iy),gd3(ix,iy));
  holdval=-surface_tension_level*curv*delta_function/(2.0*h*denval);
  px[ix][iy]+=holdval*(gd1(ix+1,iy)-gd1(ix-1,iy));
  py[ix][iy]+=holdval*(gd1(ix,iy+1)-gd1(ix,iy-1));
}   /* tensionforce */

/* compute contribution to surface tension using level set identified
   dnum */
double tangent_value_main(ix,iy,dnum,tension_constant)
int ix,iy,dnum;
double tension_constant;
{
double c1;
double delta_function,distance_dd,dx,dy,kx,ky,hx,hy,out_dist;
double hvij[2][2],kij[2][2];
int ix1,iy1;

  delta_function=(compute_delta(gdt(dnum,ix,iy))+
                  compute_delta(gdt(dnum,ix,iy+1))+
                  compute_delta(gdt(dnum,ix+1,iy+1))+
                  compute_delta(gdt(dnum,ix+1,iy)) )/4.0;
  if (delta_function>0.0)
    {
     for (ix1=0;ix1<=1;ix1++)
      for (iy1=0;iy1<=1;iy1++)
       {
        hvij[ix1][iy1]=compute_heaviside(gdt(dnum,ix+ix1,iy+iy1));
        if (stokes_breaking_flag)
         kij[ix1][iy1]=compute_surface_force(ix+ix1);
        else
         kij[ix1][iy1]=compute_curvature(ix+ix1,iy+iy1,
                         &distance_dd,&dx,&dy,&out_dist,dnum);
/*        printf("curvature ix,iy,ix1,iy1 %10.5f%5d%5d%5d%5d\n",
          kij[ix1][iy1],ix,iy,ix1,iy1);*/
        if (distance_dd>almost_same)
         {
          if ((reinit_level)&&(! hou_formulation))
           kij[ix1][iy1] /= distance_dd;
         }
        else
         kij[ix1][iy1]=0.0;
       }

     hx=(hvij[1][1]+hvij[1][0]-hvij[0][1]-hvij[0][0])/(2.0*h);
     hy=(hvij[1][1]+hvij[0][1]-hvij[1][0]-hvij[0][0])/(2.0*h);
           
     kx=(kij[1][1]+kij[1][0]-kij[0][1]-kij[0][0])/(2.0*h);
     ky=(kij[1][1]+kij[0][1]-kij[1][0]-kij[0][0])/(2.0*h);

     c1=ky*hx-kx*hy;
/*
     printf("c1,k00,h,ix,iy%10.5f%10.5f%10.5f%5d%5d\n",
        c1,kij[0][0],hvij[0][0],ix,iy);
*/
     c1*=tension_constant;
     return c1;
    }   /* delta_function>0.0 */
  else
    return 0.0;

}   /* tangent_value_main */

/* compute tension forces between fluids represented by dnum1 and dnum2 */
double tangent_value(ix,iy,dnum1,dnum2,tension_constant)
int ix,iy,dnum1,dnum2;
double tension_constant;
{
double c1,c2;

  c1 = tangent_value_main(ix,iy,dnum1,tension_constant);
  if (dnum1==dnum2)
   return c1;
  else
   {
    c2 = tangent_value_main(ix,iy,dnum2,tension_constant);
/*
if ((dnum1==1)&&(dnum2==3) && (
    ((fabs(c1)>almost_same)||(fabs(c2)>almost_same)) ))
printf("c1,c2,d1,d2,const%15.13f%15.13f%5d%5d%15.13f\n",
  c1,c2,dnum1,dnum2,tension_constant);
*/
    return (c1+c2)/2.0;
   }
}   /* tangent_value */

double check_interface(ix,iy,dmida,dmidb,dmidc,a,b,surf,cutoff) 
int ix,iy,a,b;
double dmida,dmidb,dmidc,surf,cutoff;
{
  if ((fabs(dmida)<cutoff) && (fabs(dmidb)<cutoff) && 
      (dmida*dmidb<=0.0) && (dmidc>max(dmida,dmidb)))
   {  
    if (feed_back_dnum!=-1)  
     return tangent_value(ix,iy,a,b,surf);  
    else if (b==2)  
     return tangent_value(ix,iy,a,a,surf);
    else
     return 0.0;
   }
  else
   return 0.0;
}    /* check_interface */

double get_hx(ix)
int ix;
{
  if ((ix<ixmin)&&(!use_symmetric)&&(!use_axisymmetric))
   return farx+h;
  else if (ix>=ixmax)
   return farx+h;
  else
   return h;
}

double get_hy(iy)
int iy;
{
  if (iy<iymin)
   return fary+h;
  else if (iy>=iymax)
   return fary+h;
  else
   return h;
}

/* initialize aij,bij 
   compute new value of st which is used to compute new values of ut and vt.
   px,py also updated here.
*/
void projection_fdiff(cfl_project)
double *cfl_project;
{
int ix,iy,ix1,iy1;
int iterations;
double total_tension,curvature_term,dmid1,dmid2,dmid3,
  dz,dr,divfree_ut,divfree_vt,cutoff,d1s,d2s,d3s,d4s,
  hx0,hx1,hy0,hy1,hx2,hy2,hxx0,hxx1,hyy0,hyy1;
logical drop_not_burst,symmetry_flag;
double store_density[2][2],store_beta[3][3],store_sigma[2][2],
  store_alpha[3][3],rval,rplus,rminus;

#define dist_midpoint(dnum) ((gdt(dnum,ix,iy)+gdt(dnum,ix+1,iy)+  \
                              gdt(dnum,ix,iy+1)+gdt(dnum,ix+1,iy+1))/4.0)

  drop_not_burst=((three_phase())&&(fluid_select==160)&&
   (feed_back_dnum==2)&&(tension_oil_air!=0.0)&&(tension_oil_water!=0.0));


  if (solid_body() )
   {
    for_all_points
     {
      st[ix][iy]=0.0;
      px[ix][iy]=0.0;
      py[ix][iy]=0.0;
      ut[ix][iy]=0.0;
      vt[ix][iy]=0.0;
     }

    *cfl_project=0.0;
    project_used=TRUE;
    return;
   }


  for_all_points
   if (point2_ok) /* include points outside the domain just in case periodic */
    {
     total_tension = surface_tension_level;
     if ((three_phase())&&(feed_back()) )
       total_tension += tension_oil_air + tension_oil_water;

     curvature_term = 0.0;
     if (total_tension != 0.0)
      {
       dmid1=dist_midpoint(1);
       cutoff = density_spread*h;
       if ((three_phase())&&(feed_back())&&(!drop_not_burst) )
        {
         dmid2=dist_midpoint(2);
         dmid3=dist_midpoint(3);

         curvature_term+=check_interface(ix,iy,dmid1,dmid2,dmid3,
           1,2,surface_tension_level,cutoff);
         curvature_term+=check_interface(ix,iy,dmid1,dmid3,dmid2,
           1,3,tension_oil_air,cutoff);
         curvature_term+=check_interface(ix,iy,dmid3,dmid2,dmid1,
           3,2,tension_oil_water,cutoff);
        }  /* three phase */
       else
        {  /* two phase */
         if (fabs(dmid1)<cutoff)
          curvature_term = tangent_value(ix,iy,1,1,surface_tension_level);
         if (drop_not_burst)
          {
           dmid3=dist_midpoint(3);
           if (fabs(dmid3)<cutoff)
            curvature_term+=tangent_value(ix,iy,3,3,surface_tension_level);
          }
        }
      }
     else 
      curvature_term = 0.0;  /* total_tension=0.0 */

#define get_d_val(incx,incy) cvt_density(gd1(ix+incx,iy+incy),   \
    gd2(ix+incx,iy+incy),gd3(ix+incx,iy+incy))
#define get_r_val(incx) (gr(ix+incx))

#define d_average(x1,y1,x2,y2) cvt_density(      \
     (gd1(ix+x1,iy+y1)+gd1(ix+x2,iy+y2))/2.0,    \
     (gd2(ix+x1,iy+y1)+gd2(ix+x2,iy+y2))/2.0,    \
     (gd3(ix+x1,iy+y1)+gd3(ix+x2,iy+y2))/2.0)
#define r_avg(x1,x2) ( (gr(ix+x1)+gr(ix+x2))/2.0 )

#define coeff1a(x1,x2) (get_d_val(x1,x2)/(4*get_r_val(x1)) )
#define coeff1(x1,x2) (get_d_val(x1,x2)/(2*sqr(h)*get_r_val(x1)) )
#define coeff2(x1,x2,x3,x4,x5,x6) (d_average(x1,x2,x3,x4)/(sqr(h)*r_avg(x5,x6)) )

#define force_u(incx,incy) (get_d_val(incx,incy)*gut(ix+incx,iy+incy))
#define force_v(incx,incy) (get_d_val(incx,incy)*gvt(ix+incx,iy+incy))

     a00[ix][iy]=a01[ix][iy]=a02[ix][iy]=0.0;
     a10[ix][iy]=a11[ix][iy]=a12[ix][iy]=0.0;
     a20[ix][iy]=a21[ix][iy]=a22[ix][iy]=0.0;

     hx0=hx1=hy0=hy1=h;

     if (mx_type==FREE_BDRY)
      {
       hx0=get_hx(ix-1);
       hx1=get_hx(ix);
      }
     if (my_type==FREE_BDRY)
      {
       hy0=get_hy(iy-1);
       hy1=get_hy(iy);
      }

     hx2=hy2=h;

     hxx0=hx2*hx0;hyy0=hy2*hy0;
     hxx1=hx2*hx1;hyy1=hy2*hy1;

     for (ix1=0;ix1<2;ix1++)
      for (iy1=0;iy1<2;iy1++)
       {
        store_sigma[ix1][iy1]=store_density[ix1][iy1]=
           get_d_val(ix1,iy1)/get_r_val(ix1);
       }
     for (ix1=0;ix1<3;ix1++)
      for (iy1=0;iy1<3;iy1++)
       {
        store_alpha[ix1][iy1]=store_beta[ix1][iy1]=-1.0/4.0;
       }

     if (use_corner_difference)
      {
       a00[ix][iy]=d1s=store_density[0][0]*store_beta[0][0]/hyy0+
         store_sigma[0][0]*store_alpha[0][0]/hxx0;
       a20[ix][iy]=d2s=store_density[1][0]*store_beta[2][0]/hyy0+
         store_sigma[1][0]*store_alpha[2][0]/hxx1;
       a02[ix][iy]=d3s=store_density[0][1]*store_beta[0][2]/hyy1+
         store_sigma[0][1]*store_alpha[0][2]/hxx0;
       a22[ix][iy]=d4s=store_density[1][1]*store_beta[2][2]/hyy1+
         store_sigma[1][1]*store_alpha[2][2]/hxx1;
       a01[ix][iy]=
        -store_beta[0][1]*(store_density[0][1]/hyy1+store_density[0][0]/hyy0)+
         store_alpha[0][1]*(store_sigma[0][1]/hxx0+store_sigma[0][0]/hxx0);
       a21[ix][iy]=
        -store_beta[2][1]*(store_density[1][1]/hyy1+store_density[1][0]/hyy0)+
         store_alpha[2][1]*(store_sigma[1][1]/hxx1+store_sigma[1][0]/hxx1);
       a10[ix][iy]=
         store_beta[1][0]*(store_density[1][0]/hyy0+store_density[0][0]/hyy0)-
         store_alpha[1][0]*(store_sigma[1][0]/hxx1+store_sigma[0][0]/hxx0);
       a12[ix][iy]=
         store_beta[1][2]*(store_density[1][1]/hyy1+store_density[0][1]/hyy1)-
         store_alpha[1][2]*(store_sigma[1][1]/hxx1+store_sigma[0][1]/hxx0);

       a11[ix][iy]=
         -store_beta[1][1]*(store_density[1][1]/hyy1+store_density[1][0]/hyy0+
            store_density[0][1]/hyy1+store_density[0][0]/hyy0)-
          store_alpha[1][1]*(store_sigma[1][1]/hxx1+store_sigma[1][0]/hxx1+
            store_sigma[0][1]/hxx0+store_sigma[0][0]/hxx0);
      }
     else
      {
       
       a01[ix][iy]=d1s=-coeff2(0,0,0,1,0,0);
       a10[ix][iy]=d2s=-coeff2(0,0,1,0,0,1);
       a12[ix][iy]=d3s=-coeff2(0,1,1,1,0,1);
       a21[ix][iy]=d4s=-coeff2(1,0,1,1,1,1);
       a11[ix][iy]=-(d1s+d2s+d3s+d4s);
      }

     dz=(force_u(0,1)+force_u(1,1)-force_u(0,0)-force_u(1,0))/(2.0*h);
     dr=(force_v(1,0)+force_v(1,1)-force_v(0,0)-force_v(0,1))/(2.0*h);

     bij[ix][iy]=(dr-dz)+curvature_term;
     
     if (project_used)
      st[ix][iy]=st_hld[ix][iy];
     else
      st[ix][iy]=0.0;
    }    /* point2_ok - looping all points on the STAGGERED GRID (st grid) */

#define pprintf(ss,ix,iy) symmetry_flag=FALSE
/*#define pprintf(ss,ix,iy) printf(ss,ix,iy)*/

  symmetry_flag=TRUE;
if (use_debug)
  printf("checking symmetry... \n");
  for_all_points
   if (point1_ok)
    {
     if (a12[ix][iy] != a10[ix][iy+1]) pprintf("a12 ix,iy%5d%5d\n",ix,iy);
     if (a10[ix][iy] != a12[ix][iy-1]) pprintf("a10 ix,iy%5d%5d\n",ix,iy);
     if (a22[ix][iy] != a00[ix+1][iy+1]) pprintf("a22 ix,iy%5d%5d\n",ix,iy);
     if (a00[ix][iy] != a22[ix-1][iy-1]) pprintf("a00 ix,iy%5d%5d\n",ix,iy);
     if (a02[ix][iy] != a20[ix-1][iy+1]) pprintf("a02 ix,iy%5d%5d\n",ix,iy);
     if (a20[ix][iy] != a02[ix+1][iy-1]) pprintf("a20 ix,iy%5d%5d\n",ix,iy);
     if (a01[ix][iy] != a21[ix-1][iy]) pprintf("a01 ix,iy%5d%5d\n",ix,iy);
     if (a21[ix][iy] != a01[ix+1][iy]) pprintf("a21 ix,iy%5d%5d\n",ix,iy);
    }
  if (! symmetry_flag)
   printf("bad symmetry\n");
if (use_debug)
  printf("symmetry has been checked \n");

  invert_matrix(tolerance_matrix,&iterations,st,use_corner_difference);

  *cfl_project=0.0;
  
  /* use gst3/gst2 for projection - s_t has no additional
     inflow. use gst1 for advection "s" has inflow */

  for_all_points
   if (! point1_ok)
   {
    st[ix][iy]=gst3(st,ix,iy);
   }

  for_domain
    {
     hx0=hy0=h;

     rval=1.0;
     if (use_axisymmetric)
      rval=x[ix];

     /* use gst3/gst2 for projection - s_t has no additional
        inflow. use gst1 for advection "s" has inflow */

     divfree_ut=
      (st[ix][iy]-st[ix][iy-1]+st[ix-1][iy]-st[ix-1][iy-1])/
      (2.0*hy0*rval);
     divfree_vt=
      (st[ix-1][iy-1]+st[ix-1][iy]-st[ix][iy-1]-st[ix][iy])/
      (2.0*hx0*rval);

     /* initialize (grad P)/density */
     px[ix][iy]=ut[ix][iy]-divfree_ut;
     py[ix][iy]=vt[ix][iy]-divfree_vt;
     tensionforce(ix,iy); 
     /* add surface tension effect to pressure grad/rho */

     ut[ix][iy]=divfree_ut;
     vt[ix][iy]=divfree_vt;

     *cfl_project=max(*cfl_project,
           sqrt(sqr(divfree_ut)+sqr(divfree_vt)) );

/*
printf("ut,vt,ix,iy,px,py,d%10.5f%10.5f%5d%5d%10.5f%10.5f%10.5f\n",
 ut[ix][iy],vt[ix][iy],ix,iy,px[ix][iy],py[ix][iy],d1[ix][iy]);
*/
    }
  for_all_points
   if (! point_ok)
    {
     px[ix][iy]=guu(px,ix,iy);
     py[ix][iy]=gvv(py,ix,iy);
    }
  for_all_points
   {
    px[ix][iy]*=cvt_density(gd1(ix,iy),gd2(ix,iy),gd3(ix,iy));
    py[ix][iy]*=cvt_density(gd1(ix,iy),gd2(ix,iy),gd3(ix,iy));
   }
  for (ix=0;ix<=mx;ix=ix+2)
  {
   iy=0;
   if (ix==0) 
    pp[ix][iy]=0.0;
   else  
    pp[ix][iy]=pp[ix-1][iy]+h*(px[ix][iy]+px[ix][iy+1])/2.0;
   if (ix<mx)
    pp[ix+1][iy]=pp[ix][iy]+h*(px[ix+1][iy]+px[ix+1][iy+1])/2.0;

   for (iy=1;iy<=my;iy++)
    {
     pp[ix][iy]=pp[ix][iy-1]+h*(py[ix][iy]+py[ix+1][iy])/2.0;
     if (ix<mx)
      pp[ix+1][iy]=pp[ix][iy]+h*(px[ix+1][iy]+px[ix+1][iy+1])/2.0;
    }
  }   /* end ix loop */

  for_all_points
   {
    px[ix][iy]/=cvt_density(gd1(ix,iy),gd2(ix,iy),gd3(ix,iy));
    py[ix][iy]/=cvt_density(gd1(ix,iy),gd2(ix,iy),gd3(ix,iy));
   }

  project_used=TRUE;
}   /* projection_fdiff */
